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Abstract 



We present an approach for the joint segmentation and grouping of similar com- 
ponents in anisotropic 3D image data and use it to segment neural tissue in serial 
sections electron microscopy (EM) images. We first construct a nested set of neu- 
ron segmentation hypotheses for each slice. A conditional random field (CRF) 
then allows us to evaluate both the compatibility of a specific segmentation and 
a specific inter-slice assignment of neuron candidates with the underlying obser- 
vations. The model is solved optimally for an entire image stack simultaneously 
using integer linear programming (ILP), which yields the maximum a posteriori 
solution in amortized linear time in the number of slices. We evaluate the perfor- 
mance of our approach on an annotated sample of the Drosophila larva neuropil 
and show that the consideration of different segmentation hypotheses in each slice 
leads to a significant improvement in the segmentation and assignment accuracy. 



1 Introduction 



Electron microscopy (EM) remains the only imaging technique with sufficient resolution for the 
elucidation of synaptic contacts between neurons (3). The acquisition of large volumes of brain 
circuitry is now possible with recent advances in automated imaging for EM. The next bottleneck in 
the reconstruction of neural circuits is the accurate reconstruction of 3D neural arbors from stacks of 
EM images |5 ]. Despite several efforts at automating the analysis of these highly stereotyped images, 
further improvements in overall accuracy are needed before automated methods can eliminate the 
tedious work currently needed to annotate large EM image datasets. 

Current attempts at automatic labelling of anisotropic neural tissue in 3D-imaged serial sections can 
broadly be divided in two categories, (i) Binary segmentation based approaches try to identify the 
outlines of neurons within single slices (7}|9), and the results are used to establish geometrically con- 
sistent assignments of connected components that belong to one neuron. These approaches suffer 
from a high sensitivity to small errors: a missing piece of membrane alters the topological proper- 
ties of the result, (ii) Over- segmentation based approaches merge small image regions within and 
between slices (14][T5). The resulting optimization problem is solved approximately. 
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Figure 1: The proposed processing pipeline. 



Our approach belongs to the first category. The improvements over current work are: (i) Our frame- 
work allows for rivaling concurrent segmentation hypotheses. This reduces the likelihood of missing 
crucial parts of the segmentation, (ii) All possible continuations of neural processes across slices are 
considered jointly with the segmentation hypotheses. Thus, higher-order relations between slices 
influence the segmentation, (iii) A globally consistent and optimal solution is found in amortized 
linear time in the number of slices. This is a prerequisite for the processing of large volumes. 

An overview of our approach can be seen in Fig. [T] From the original images, per-pixel local 
features are extracted. Using several segmentations based on predictions of a pre-trained random 
forest classifier on the local features, segmentation hypotheses are extracted for each slice. These 
hypotheses are connected components of pixels (see Sec. [2]). Given the hypotheses, two tasks have 
to be performed. Firstly, a subset of non-overlapping components has to be found to make up the 
segmentation of each slice. Secondly, assignments have to be established between the selected 
components of two subsequent slices. These assignments identify components that belong to the 
same neuron. We show how both tasks can be solved jointly by the introduction of binary assignment 
variables (see Sec. [5]). A consistent and optimal segmentation and assignment can be found by MAP- 
inference in a conditional random field (CRF). Results presented in Sec. [4] show that our approach 
leads to a significant improvement in the segmentation and assignment accuracy. 



2 Generating Segmentation Hypotheses 



To extract segmentation hypotheses for each slice, we use a parameterized model to solve a binary 
image segmentation task. Given a field of per-pixel local features x, the task is to maximize the 
conditional probability over a binary segmentation y: 



p(y\x) 
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Here, ft c M? is the image domain and D(xi,yi) = logp(xi\yi) is the log-likelihood of pixel 
i e ft belonging to fore- or background, as given by a pre-trained random forest classifier. The set 
AT C ft x ft contains all pairs of 8-connected neighboring pixels. Z(x) is the partition function. The 
second term in the exponent ensures smoothness of the segmentation, while favoring label changes 
in image regions of strong spatial gradients fT0| : 
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Here, g(xi, Xj) measures the difference in gray- leveljQ ||.|| denotes the length of a vector, and 5 is 
the Kronecker-delta. a is a parameter of the smoothness term. 

The last term in the exponent of ([T]) is a prior on the expected number of pixels assigned to the 
label 'neuron'. The scalars A^, As, and A at are parameters of the probability distribution. Since 
the exponent in ([T]) is submodular, the inference task can be seen as a parametric max-flow problem 
which can be solved efficiently (2). However, finding an optimal set of parameters is a non-trivial 
task and one cannot expect a fixed set of parameters to perform well on all images (TTJ. We go further 
and claim that a fixed set of parameters cannot even be expected to perform well on all areas of one 
image. Therefore, we enumerate several local segmentation hypotheses by variation of A at (Fig. [2]). 



! For simplicity, we assume that the gray level is part of the local feature vector x% 
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Figure 3: Visualisation of the segmentation hypotheses extraction. For different values of the prior 
parameter A at (shades of blue) connected components of the segmentation are found (left side). The 
subset relation of these connected components define the component tree (right side). 



This can be done efficiently by several warm-started graph-cuts fTQ| , from which we obtain a series 
of segmentations. Each segmentation consists of a set of connected components C % C ft that are 
labelled as neurons. Each of these components is considered as one segmentation hypothesis. As A at 
decreases, these components can grow and merge, thus establishing a tree-shaped subset hierarchy, 
the so-called component tree |6 ] (Fig. [3]). All components which do not meet a criterion for stability 
over the range of values of A at are removed (12). Let C z denote the set of all remaining components 
of slice z. Any consistent subset S z C C z of these hypotheses yields a valid segmentation of this 
slice. A subset is consistent if none of the containing components overlap, i.e., C l D C J = 
for all C l ,C J G S z with i ^ j. In the following, we show how we find the optimal consistent 
segmentation by considering the assignments of hypotheses between pairs of adjacent slices. 



3 Assignment Model 

For each possible assignment of a segmentation hypothesis in one slice to a hypothesis in the pre- 
vious or next slice, we introduce one binary assignment variable. This variable is set to 1 if the 
involved hypotheses and their mutual assignment are accepted. 

Let m be the number of all possible assignments. A binary vector a G {0, l} m of assignment vari- 
ables is created similarly to the method proposed in 1 13]. Each possible continuation of a hypothesis 
C % in slice z to a hypothesis C J in slice z + 1 is represented by a variable a*^ J . A split of C l in slice 
z to C- 7 and C k in slice z + 1 is represented as a l ^ ,k . Similarly, each possible merge is encoded 
as a l ^^ k . Appearances and disappearances of neurons are encoded as assignments to a special end 
node E, i.e., for each hypothesis C % we introduce two variables a l ^ E and a E ^ % . For the possible 
assignments, only components within a threshold distance to each other are considered. Thus, the 
number of assignment variables is linear in the number of segmentation hypotheses. See Fig. [4] for 
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Figure 4: Examples of the four outgoing assignment categories for a single segmentation hypoth- 
esis (orange): A continuation (red) is modelled for each hypothesis in the next slice that is within 
a threshold distance. Possible splits (green) and merges (blue) are enumerated for neighboring hy- 
potheses in the respective slices. Again, the possible sources or targets of splits and merges are 
segmentation hypotheses within a threshold distance. The disappearance of a neuron is represented 
by a single assignment for each candidate (black). 



examples of assignments of a single segmentation hypothesis. We find the optimal segmentation 
and assignment via MAP-inference on the following CRF: 



p(a|c) = zk exp 



-c T a- £ff(P)- Y, 

Per i<i<n 



(3) 



Here, V denotes the set of all complete paths of the hypotheses trees. The set of all incoming 
and outgoing assignment variables for a hypothesis C % is denoted by a^* and a 2 ^, respectively. 
Incoming assignment variables of a hypothesis C l are all components of a that have C l on the right 
side of the superscript. Outgoing assignment variables are defined analogously. 

The first term in the exponent accounts for unary potentials determining the costs for each assign- 
ment. For that, a vector c G M m is constructed, which is congruent to a. The costs for one-to-one 
assignments are modelled as: 

= L L(C ij ) + P \\& - C^ll 2 + 5 1 C* C j \ 2 . (4) 

Here, we write C lJ as a shorthand for C l U C J . The expression C l denotes the mean pixel position 
of the candidate C % ; C % C 3 is the mean-corrected symmetric set difference of candidates C % 
and C 3 \ ||.|| denotes the Euclidean distance; and |.| the cardinality of a set. The term L(C) is 
proportional to the probability of assigning all pixels of C to a neuron, i.e., 

L(C) = ^2(D(x j ,0)-D(x j ,l)) + Sj, k (y j: y k ,x). (5) 

jec, k^c 

In a similar way, we define the costs for splits: 

c i^j,k = 0LL ( C v k ) + e BP \\c l - c^\\ 2 + 6 B s\c l e c 3k \ 2 . (6) 

The merge cases are defined analogously. Costs for the appearance or disappearance of a neuron 
depend on the data term and the size of the component: 

C i-+E = c E^i = 0lL ^ + e E \&\ 2 . (7) 

The scalars Ol, Op, Os, Obp, Obs, and 6e are parameters of the model. 

The second and third terms in ^ are higher-order potentials that ensure the consistency of the so- 
lution. H(P) ensures that at most one of all incoming assignments of competing segmentation 
hypotheses is selected. Competing hypotheses are components that share some pixels, i.e., all com- 



ponents along one path in the hypotheses tree (see Fig. 5(a)). This hypothesis constraint ensures that 
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(a) hypothesis consistency (b) explanation consistency 



Figure 5: Visualization of the two types of consistency constraints. The hypothesis consistency (a) 
ensures that no pixel is assigned to more than one segmentation hypothesis: For each path of the 
component tree (blue), the sum of all incoming assignment variables (gray) has to at most one. The 



explanation consistency (b) ensures a continuous sequence of assignments: For each segmentation 
hypothesis (orange), the sum of all incoming assignment variables (from the previous image) has 
to be equal to the sum of all outgoing assignment variables (to the next image). The incoming or 
outgoing assignment variables for a component are all assignment variables that have the component 
as target or source, respectively. 



each pixel is explained by at most one hypothesis and that at most one of the possible assignments 
is selected for each hypothesis. 



H(P) 



if £ i€ pEa€a-««<l 

oo else 



(8) 



E(8l~^ % , a 1 ^) ensures that if an incoming assignment was selected for a segmentation hypothesis, 
an outgoing assignment is select ed as well. This explanation constraint guarantees a consistent 
sequence of assignments (see Fig. [5(b)). 
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oo else 



(9) 



Since the potentials H and E impose hard constraints and the remaining potentials are linear in a, 
the MAP solution to ^ can be found by solving the following linear program. 

minc T a s.t. a - X VPgP (10) 

a ~ a ' = l<i<n (11) 

Unfortunately, the constraint matrix of this linear program is not totally unimodular. Therefore, we 
have to enforce the integrality of the solution explicitly. The resulting optimization problem is an 
instance of an integer linear program (ILP), which we solve using the IBM CPLEX solver (TJ. 



4 Results 



We evaluated the performance of our approach on an annotated sample of Drosophila larva neu- 
ral tissue fit. This publicly available data set consists of 30 serial sections (50 nm), imaged with 
transmission electron microscopy at a resolution of 4x4x50 nm/pixel. The image volume contains 
a 2x2x1.5 micron cube of neuropil tissue. The dataset includes labels of cellular membranes, cyto- 
plasms and mitochondria of all 170 neural processes contained in the dataset. 

To measure the accuracy of our method, we use the edit distance between the result and the ground- 
truth, i.e., the number of splits and merges a human operator would need to perform to restore the 
ground-truth |T5| . In contrast to the measure proposed in (5|, we count every false merge, even if the 
same objects are involved. Each missed neuron segment in a slice is counted as one merge error and 
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Figure 6: Comparison of the accuracy of our approach with different segmentation hypotheses (black 
square) against a series of experiments without segmentation hypotheses (colored circles). The con- 
sideration of different segmentations is superior to all possible single segmentations. The single 
segmentation experiments cover the whole range from under- segmentation (blue) to over segmenta- 
tion (red). 



each falsely introduced neuron segment is counted as one split error. Between the slices, we count 
each missed assignment as a split error and each falsely introduced assignment as a merge error. 

To evaluate the performance of our approach in comparison to approaches that do not allow local 
variations in the segmentation parameters (which is the case for all existing approaches that we 
know of), we carried out a series of experiments with fixed A at (see Eq.[T]). In particular, we took 
34 equidistant samples of A at in an interval reaching from obvious over- to under- segmentation of 
the images. Note that the clamping of the parameter A at corresponds to not considering conflicting 
hypotheses: all the component trees have a depth of one. All other parameters of the pipeline do not 
change between the experiments, and have been found via a grid-search in the parameter space. We 
compare the results of this series with a single run of our proposed method that does allow for local 
variations of A at. In Fig. [6] we show the edit distance for each A n in the series as well as for our 
approach normalized by the number of neurons in the dataset. Our data demonstrates that the single 
run with local variation of A at is superior to all the cases in which A at was fixed. 

In Fig. [7] we show a representative segmentation example of five subsequent slices using our ap- 
proach. In the same figure we give the inference time for the processing of the test dataset with 
different number of slices. The results indicate that the solution to the optimization task can be 
found in amortized linear time. 



5 Conclusion 

We presented a novel approach for the joint segmentation and assignment of similar regions in 
anisotropic 3D image data with an amortized linear running time in the number of slices. On 
an annotated sample of neural tissue we could show that the consideration of different segmen- 
tation hypotheses significantly improves the accuracy of the segmentation and grouping of neural 
processes. The probabilistic formulation of our approach is particularly useful in semi-interactive 
environments, where human-reconstructed neurons can help resolving ambiguity. Every manually 
traced neural process can be used to impose strong priors on the assignment model by decreasing 
the costs for the respective assignments. 
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Figure 7: Sample of the segmentation result of five subsequent slices (a) and inference time (b). 
Examples of each error type are highlighted. Colors indicate the assignment of components between 
the slices. 
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Currently, the assignment costs are computed from simple statistics of the involved segmentation 
hypotheses. Improvements might be obtained by using more sophisticated features (8). This will be 
the focus of further research. 

6 Acknowledgements 

This work was funded by the Swiss National Science Foundation. Parts of this work have been 
completed in collaboration with the Heidelberg Collaboratory for Image Processing. 

References 

[I] IBM ILOG CPLEX Optimizer, v!2.2. |http : //www- 01 . ibm. com/sof tware/integration/| 
|opt imizat ion/ cp lex- opt imizer/| 

[2] Yuri Boykov and Vladimir Kolmogorov, An experimental comparison of min- cut/max- flow algorithms 
for energy minimization in vision, IEEE Transactions on Pattern Analysis and Machine Intelligence 26 
(2004sept.), no. 9, 1 124 -1 137. 

[3] Kevin L Briggman and Winfried Denk, Towards neural circuit reconstruction with volume electron mi- 
croscopy techniques, Current Opinion in Neurobiology 16 (2006), no. 5, 562 -570. 

[4] Albert Cardona, Stephan Saalfeld, Stephan Preibisch, Benjamin Schmid, Anchi Cheng, Jim Pulokas, 
Pavel Tomancak, and Volker Hartenstein, An integrated micro- and macroarchitectural analysis of the 
Drosophila brain by computer-assisted serial section electron microscopy, PLoS Biology 8 (2010), no. 10, 
el00050. 

[5] Viren Jain, H. Sebastian Seung, and Srinivas C. Turaga, Machines that learn to segment images: a crucial 
technology for connectomics, Current Opinion in Neurobiology (2010). 

[6] Component Trees for Image Filtering and Segmentation (1997) 

[7] Elizabeth Jurrus, Antonio R.C. Paiva, Shigeki Watanabe, James R. Anderson, Bryan W. Jones, Ross T. 
Whitaker, Erik M. Jorgensen, Robert E. Marc, and Tolga Tasdizen, Detection of neuron membranes 
in electron microscopy images using a serial neural network architecture, Medical Image Analysis 14 
(2010), no. 6,770-783. 

[8] Verena Kaynig, Thomas Fuchs, and Joachim M. Buhmann, Geometrical Consistent 3D Tracing of 
Neuronal Processes in ssTEM Data, Proceedings of the Conference on Medial Image Computing and 
Computer- Assisted Intervention, 2010, pp. 209-216. 

[9] , Neuron Geometry Extraction by Perceptual Grouping in ssTEM Images, Proceedings of the IEEE 

Conference on Computer Vision and Pattern Recognition, 2010, pp. 2902-2909. 

[10] Pushmeet Kohli and Philip Torr, Dynamic Graph Cuts and Their Applications in Computer Vision, Com- 
puter Vision, 2010, pp. 51-108. 

[II] Vladimir Kolmogorov, Yuri Boykov, and Cars ten Rother, Applications of Parametric Maxflow in Com- 
puter Vision, Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2007. 

[12] Jiri Matas, Ondrej Chum, Martin Urban, and Tomas Pajdla, Robust Wide-Baseline Stereo from Maximally 
Stable Extremal Regions, Image and Vision Computing 22 (2004), no. 10, 761 -767. 

[13] Dirk Padfield, Jens Rittscher, and Badrinath Roysam, Coupled minimum-cost flow cell tracking for high- 
throughput quantitative analysis, Medical Image Analysis In Press, Corrected Proof (2010). 

[14] Amelio Vazquez-Reina, Shai Avidan, Hanspeter Pfister, and Eric Miller, Multiple Hypothesis Video Seg- 
mentation from Superpixel Plows, Proceedings of the European Conference on Computer Vision (ECCV), 
2010. 

[15] Shiv Naga Prasad Vitaladevuni and Ronen Basri, Co-Clustering of Image Segments Using Convex Op- 
timization Applied to EM Neuronal Reconstruction, Proceedings of the IEEE Conference on Computer 
Vision and Pattern Recognition, 2010, pp. 2203 -2210. 



8 



